home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / integration / qawf.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-22  |  5.8 KB  |  282 lines
  1. /* integration/qawf.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <math.h>
  22. #include <float.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_integration.h>
  26.  
  27. #include "initialise.c"
  28. #include "append.c"
  29. #include "qelg.c"
  30.  
  31. int
  32. gsl_integration_qawf (gsl_function * f,
  33.               const double a,
  34.               const double epsabs,
  35.               const size_t limit,
  36.               gsl_integration_workspace * workspace,
  37.               gsl_integration_workspace * cycle_workspace,
  38.               gsl_integration_qawo_table * wf,
  39.               double *result, double *abserr)
  40. {
  41.   double area, errsum;
  42.   double res_ext, err_ext;
  43.   double correc, total_error = 0.0, truncation_error;
  44.  
  45.   size_t ktmin = 0;
  46.   size_t iteration = 0;
  47.  
  48.   struct extrapolation_table table;
  49.  
  50.   double cycle;
  51.   double omega = wf->omega;
  52.  
  53.   const double p = 0.9;
  54.   double factor = 1;
  55.   double initial_eps, eps;
  56.   int error_type = 0;
  57.  
  58.   /* Initialize results */
  59.  
  60.   initialise (workspace, a, a);
  61.  
  62.   *result = 0;
  63.   *abserr = 0;
  64.  
  65.   if (limit > workspace->limit)
  66.     {
  67.       GSL_ERROR ("iteration limit exceeds available workspace", GSL_EINVAL) ;
  68.     }
  69.  
  70.   /* Test on accuracy */
  71.  
  72.   if (epsabs <= 0)
  73.     {
  74.       GSL_ERROR ("absolute tolerance epsabs must be positive", GSL_EBADTOL) ;
  75.     }
  76.  
  77.   if (omega == 0.0)
  78.     {
  79.       if (wf->sine == GSL_INTEG_SINE)
  80.     {
  81.       /* The function sin(w x) f(x) is always zero for w = 0 */
  82.  
  83.       *result = 0;
  84.       *abserr = 0;
  85.  
  86.       return GSL_SUCCESS;
  87.     }
  88.       else
  89.     {
  90.       /* The function cos(w x) f(x) is always f(x) for w = 0 */
  91.  
  92.       int status = gsl_integration_qagiu (f, a, epsabs, 0.0,
  93.                           cycle_workspace->limit,
  94.                           cycle_workspace,
  95.                           result, abserr);
  96.       return status;
  97.     }
  98.     }
  99.  
  100.   if (epsabs > GSL_DBL_MIN / (1 - p))
  101.     {
  102.       eps = epsabs * (1 - p);
  103.     }
  104.   else
  105.     {
  106.       eps = epsabs;
  107.     }
  108.  
  109.   initial_eps = eps;
  110.  
  111.   area = 0;
  112.   errsum = 0;
  113.  
  114.   res_ext = 0;
  115.   err_ext = GSL_DBL_MAX;
  116.   correc = 0;
  117.  
  118.   cycle = (2 * floor (fabs (omega)) + 1) * M_PI / fabs (omega);
  119.  
  120.   gsl_integration_qawo_table_set_length (wf, cycle);
  121.  
  122.   initialise_table (&table);
  123.  
  124.   for (iteration = 0; iteration < limit; iteration++)
  125.     {
  126.       double area1, error1, reseps, erreps;
  127.  
  128.       double a1 = a + iteration * cycle;
  129.       double b1 = a1 + cycle;
  130.  
  131.       double epsabs1 = eps * factor;
  132.  
  133.       int status = gsl_integration_qawo (f, a1, epsabs1, 0.0, limit,
  134.                      cycle_workspace, wf,
  135.                      &area1, &error1);
  136.  
  137.       append_interval (workspace, a1, b1, area1, error1);
  138.  
  139.       factor *= p;
  140.  
  141.       area = area + area1;
  142.       errsum = errsum + error1;
  143.  
  144.       /* estimate the truncation error as 50 times the final term */
  145.  
  146.       truncation_error = 50 * fabs (area1);
  147.  
  148.       total_error = errsum + truncation_error;
  149.  
  150.       if (total_error < epsabs && iteration > 4)
  151.     {
  152.       goto compute_result;
  153.     }
  154.  
  155.       if (error1 > correc)
  156.     {
  157.       correc = error1;
  158.     }
  159.  
  160.       if (status)
  161.     {
  162.       eps = GSL_MAX_DBL (initial_eps, correc * (1.0 - p));
  163.     }
  164.  
  165.       if (status && total_error < 10 * correc && iteration > 3)
  166.     {
  167.       goto compute_result;
  168.     }
  169.  
  170.       append_table (&table, area);
  171.  
  172.       if (table.n < 2)
  173.     {
  174.       continue;
  175.     }
  176.  
  177.       qelg (&table, &reseps, &erreps);
  178.  
  179.       ktmin++;
  180.  
  181.       if (ktmin >= 15 && err_ext < 0.001 * total_error)
  182.     {
  183.       error_type = 4;
  184.     }
  185.  
  186.       if (erreps < err_ext)
  187.     {
  188.       ktmin = 0;
  189.       err_ext = erreps;
  190.       res_ext = reseps;
  191.  
  192.       if (err_ext + 10 * correc <= epsabs)
  193.         break;
  194.       if (err_ext <= epsabs && 10 * correc >= epsabs)
  195.         break;
  196.     }
  197.  
  198.     }
  199.  
  200.   if (iteration == limit)
  201.     error_type = 1;
  202.  
  203.   if (err_ext == GSL_DBL_MAX)
  204.     goto compute_result;
  205.  
  206.   err_ext = err_ext + 10 * correc;
  207.  
  208.   *result = res_ext;
  209.   *abserr = err_ext;
  210.  
  211.   if (error_type == 0)
  212.     {
  213.       return GSL_SUCCESS ;
  214.     }
  215.  
  216.   if (res_ext != 0.0 && area != 0.0)
  217.     {
  218.       if (err_ext / fabs (res_ext) > errsum / fabs (area))
  219.     goto compute_result;
  220.     }
  221.   else if (err_ext > errsum)
  222.     {
  223.       goto compute_result;
  224.     }
  225.   else if (area == 0.0)
  226.     {
  227.       goto return_error;
  228.     }
  229.  
  230.   if (error_type == 4)
  231.     {
  232.       err_ext = err_ext + truncation_error;
  233.     }
  234.  
  235.   goto return_error;
  236.  
  237. compute_result:
  238.  
  239.   *result = area;
  240.   *abserr = total_error;
  241.  
  242. return_error:
  243.  
  244.   if (error_type > 2)
  245.     error_type--;
  246.  
  247.   if (error_type == 0)
  248.     {
  249.       return GSL_SUCCESS;
  250.     }
  251.   else if (error_type == 1)
  252.     {
  253.       GSL_ERROR ("number of iterations was insufficient", GSL_EMAXITER);
  254.     }
  255.   else if (error_type == 2)
  256.     {
  257.       GSL_ERROR ("cannot reach tolerance because of roundoff error",
  258.          GSL_EROUND);
  259.     }
  260.   else if (error_type == 3)
  261.     {
  262.       GSL_ERROR ("bad integrand behavior found in the integration interval",
  263.          GSL_ESING);
  264.     }
  265.   else if (error_type == 4)
  266.     {
  267.       GSL_ERROR ("roundoff error detected in the extrapolation table",
  268.          GSL_EROUND);
  269.     }
  270.   else if (error_type == 5)
  271.     {
  272.       GSL_ERROR ("integral is divergent, or slowly convergent",
  273.          GSL_EDIVERGE);
  274.     }
  275.   else
  276.     {
  277.       GSL_ERROR ("could not integrate function", GSL_EFAILED);
  278.     }
  279.  
  280. }
  281.  
  282.